import numpy as np
from scipy import integrate
from scipy.optimize import brentq,fsolve, root,bisect
from sympy import diff, log
import sympy as sp
import matplotlib.pyplot as plt
import pandas as pd
# import PE_cal
import traceback

# 相互作用系数 La Ce Gd Y Nd
sys = ['Mg_Al','Mg_La','Mg_Ce','Mg_Gd','Mg_Y','Mg_Nd','Al_La','Al_Ce','Al_Gd','Al_Y','Al_Nd']

def parameters(system,T):
    try:
            
        L12,theta,M,rou=[],[],[],[]
        # theta M rou
        Mg = [(577-0.26*(T-923))/1000, 0.024305, 1740]
        Al = [(871-0.155*(T-933))/1000, 0.026982, 2700]
        Zn = [(789-0.21*(T-693))/1000, 0.06539, 7140]
        La = [(728-0.1*(T-1193))/1000, 0.138906, 6150]
        Ce = [(707-0.078*(T-1077))/1000, 0.140116, 6770]
        Gd = [(664-0.058*(T-1625))/1000, 0.15725, 7900]
        Y = [(872-0.086*(T-1793))/1000, 0.088906, 4470]
        Nd = [(687-0.0885*(T-1297))/1000, 0.14424, 7010]
        # 确保T是浮点数
        T = float(T)
        
        if system == 'Mg_Al':
            L12 = [-12000+T*8.566, 1890-3*T, 2000, 0]
            theta = [Mg[0], Al[0]]
            M = [Mg[1], Al[1]]
            rou = [Mg[2], Al[2]]
        elif system == 'Mg_La':
            L12 = [-33288.5+T*8.9095, -(41081-29.4287*T), -12051.9, 0]
            theta = [Mg[0], La[0]]
            M = [Mg[1], La[1]]
            rou = [Mg[2], La[2]]
        elif system == 'Mg_Ce':
            L12 = [-39381.2+T*16.3405, -25338.6, -15106.9, 0]
            theta = [Mg[0], Ce[0]]
            M = [Mg[1], Ce[1]]
            rou = [Mg[2], Ce[2]]
        elif system == 'Mg_Gd':
            L12 = [-36681.3+T*16.2484, -(34233.8-10.7783*T), -7352.9, 0]
            theta = [Mg[0], Gd[0]]
            M = [Mg[1], Gd[1]]
            rou = [Mg[2], Gd[2]]
        elif system == 'Mg_Y':
            L12 = [-40917.97+T*22.83, -18685.23+10.24*T, 1076.67-6.05*T, 0]
            theta = [Mg[0], Y[0]]
            M = [Mg[1], Y[1]]
            rou = [Mg[2], Y[2]]
        elif system == 'Mg_Nd':
            L12 = [-34308.7+T*10.9, 14565.6, -16479, 0]
            theta = [Mg[0], Nd[0]]
            M = [Mg[1], Nd[1]]
            rou = [Mg[2], Nd[2]]
        elif system == 'Al_La':
            L12 = [-89252, 0, 0, 0]
            theta = [Al[0], La[0]]
            M = [Al[1], La[1]]
            rou = [Al[2], La[2]]
        elif system == 'Al_Ce':
            L12 = [-141000+43.82*T, -55600+23.68*T, -1686, 0]
            theta = [Al[0], Ce[0]]
            M = [Al[1], Ce[1]]
            rou = [Al[2], Ce[2]]
        elif system == 'Al_Gd':
            L12 = [-166500+52.36*T, -23790, 18520, 0]
            theta = [Al[0], Gd[0]]
            M = [Al[1], Gd[1]]
            rou = [Al[2], Gd[2]]
        elif system == 'Al_Y':
            L12 = [-160876.36+5*T, -32000+7.56*T, 32000-6.53*T, 0]
            theta = [Al[0], Y[0]]
            M = [Al[1], Y[1]]
            rou = [Al[2], Y[2]]
        elif system == 'Al_Nd':
            L12 = [-145100+44.97*T, -28030, -6708, 0]
            theta = [Al[0], Nd[0]]
            M = [Al[1], Nd[1]]
            rou = [Al[2], Nd[2]]
        elif system == 'Mg_Zn':
            L12 = [-77729.2+T*680.523-95*T*log(T)+0.04*T**2,1890-3*T,2000,0]
            theta = [Mg[0],Zn[0]]
            M = [Mg[1],Zn[1]]
            rou = [Mg[2],Zn[2]]
        elif system == 'Zn_La':
            L12 = [-89252,0,0,0]
            theta = [Zn[0],La[0]]
            M = [Zn[1],La[1]]
            rou = [Zn[2],La[2]]
        elif system == 'Zn_Ce':
            L12 = [-141000+43.82*T,-55600+23.68*T,-1686,0]
            theta = [Zn[0],Ce[0]]
            M = [Zn[1],Ce[1]]
            rou = [Zn[2],Ce[2]]
        elif system == 'Zn_Gd':
            L12 = [-166500+52.36*T,-23790,18520,0]
            theta = [Zn[0],Gd[0]]
            M = [Zn[1],Gd[1]]
            rou = [Zn[2],Gd[2]]
        elif system == 'Zn_Y':
            L12 = [-160876.36+5*T,-32000+7.56*T,32000-6.53*T,0]
            theta = [Zn[0],Y[0]]
            M = [Zn[1],Y[1]]
            rou = [Zn[2],Y[2]]
        elif system == 'Zn_Nd':
            L12 = [-145100+44.97*T,-28030,-6708,0]
            theta = [Zn[0],Nd[0]]
            M = [Zn[1],Nd[1]]
            rou = [Zn[2],Nd[2]]
        else:
            raise ValueError(f"Unsupported system: {system}")

        # 确保所有参数都是浮点数
        L12 = [float(x) for x in L12]
        theta = [float(x) for x in theta]
        M = [float(x) for x in M]
        rou = [float(x) for x in rou]

        return L12, theta, M, rou
        
    except Exception as e:
        print(f"Error in parameters function: {str(e)}")
        print(f"Stack trace: {traceback.format_exc()}")
        raise

def butler(T,beta,l_12,theta,M,rou,x1):
    try:
        # 确保所有输入都是浮点数
        T = float(T)
        beta = float(beta)
        x1 = float(x1)
        
        # 确保列表中的元素都是浮点数，并检查长度
        l_12 = [float(x) for x in l_12]
        if len(l_12) < 4:
            l_12.extend([0.0] * (4 - len(l_12)))  # 确保l_12有4个元素
            
        theta = [float(x) for x in theta]
        M = [float(x) for x in M]
        rou = [float(x) for x in rou]
        
        # 打印调试信息
        print(f"Input parameters after conversion:")
        print(f"T: {T}, beta: {beta}, x1: {x1}")
        print(f"l_12: {l_12}")
        print(f"theta: {theta}")
        print(f"M: {M}")
        print(f"rou: {rou}")
        
        R = 8.314
        N0 = 6.02E+23
        # 二元 体系
        L = 1.091
        # 使用numpy的power函数替代pow
        A = [L * np.power(N0,1/3) * np.power(M[0]/rou[0], 2/3), 
             L * np.power(N0, 1/3) * np.power(M[1]/rou[1], 2/3)]
        
        # 打印中间计算结果
        print(f"Intermediate calculations:")
        print(f"A: {A}")
        
        # deltaG 表达式
        x2 = 1-x1
        # 使用numpy的power函数
        deltaG = x1 * x2 * (l_12[0] + l_12[1] * (x1 - x2) + 
                           l_12[2] * np.power(x1-x2, 2) + 
                           l_12[3] * np.power(x1-x2, 3))
        
        deltaG_B1 = deltaG + (1-x1)*((1-2*x1)*(l_12[0]+l_12[1]*(2*x1-1)+
                        l_12[2]*np.power(2*x1-1,2)+
                        l_12[3]*np.power(2*x1-1,3))+
                x1*(1-x1)*(2*l_12[1]+4*l_12[2]*(2*x1-1)+
                        6*l_12[3]*np.power(2*x1-1,2)))
        
        deltaG_B2 = deltaG + (1-x2)*((1-2*x2)*(l_12[0]+l_12[1]*(1-2*x2)+
                        l_12[2]*np.power(2*x2-1,2)+
                        l_12[3]*np.power(2*x2-1,3))+
                x2*(1-x2)*(-2*l_12[1]-4*l_12[2]*(1-2*x2)-
                        6*l_12[3]*np.power(1-2*x2,2)))

        # print(f"deltaG calculations:")
        # print(f"deltaG: {deltaG}")
        # print(f"deltaG_B1: {deltaG_B1}")
        # print(f"deltaG_B2: {deltaG_B2}")

        def equation(x):
            x = float(x)
            # 添加数值稳定性检查
            if x <= 0 or x >= 1:
                return float('inf')
                
            # 使用numpy的log函数，并添加数值稳定性检查
            log_x_x1 = np.log(np.maximum(x / x1, 1e-300))
            log_1x_x2 = np.log(np.maximum((1 - x) / x2, 1e-300))
            
            term1 = theta[0] + (R * T * log_x_x1 + (beta-1) * deltaG_B1) / A[0]
            term2 = theta[1] + (R * T * log_1x_x2 + (beta-1) * deltaG_B2) / A[1]
            result = term1 - term2
            
            print(f"Equation evaluation at x={x}: term1={term1}, term2={term2}, result={result}")
            return result

        # 使用更安全的边界值
        try:
            solution = bisect(equation, 1e-4, 1 - 1e-4)
            print(f"Solution found: {solution}")
        except ValueError as e:
            print(f"Error in bisect: {str(e)}")
            raise

        # 使用numpy的log函数，并添加数值稳定性检查
        log_solution_x1 = np.log(np.maximum(solution / x1, 1e-300))
        theta12 = theta[0] + (R * T * log_solution_x1 + (beta-1) * deltaG_B1) / A[0]

        # 限制输出的小数位数
        theta12 = round(theta12, 6)
        print(f"Final result: {theta12}")
        return theta12
        
    except Exception as e:
        print(f"Error in butler function: {str(e)}")
        print(f"Input parameters: T={T}, beta={beta}, l_12={l_12}, theta={theta}, M={M}, rou={rou}, x1={x1}")
        print(f"Stack trace: {traceback.format_exc()}")
        raise


def calculate_butler_3(T,beta,l_12,l_23,l_13,theta,M,rou,x1,x2,x3):

    l_31 = [l_13[0], -l_13[1], l_13[2], -l_13[3]]

    # Define sigma functions
    def sigma_E12(x):
        return x * (1 - x) * (l_12[0] + l_12[1] * (2 * x - 1) + l_12[2] * (2 * x - 1) ** 2 + l_12[3]*pow((2*x-1),3))
    def sigma_E13(x):
        return x * (1 - x) * (l_31[0] - l_31[1] * (2 * x - 1) + l_31[2] * (2 * x - 1) ** 2 - l_31[3]*pow((2*x-1),3))
    def sigma_E21(x):
        return x * (1 - x) * (l_12[0] - l_12[1] * (2 * x - 1) + l_12[2] * (2 * x - 1) ** 2 - l_12[3]*pow((2*x-1),3))
    def sigma_E23(x):
        return x * (1 - x) * (l_23[0] + l_23[1] * (2 * x - 1) + l_23[2] * (2 * x - 1) ** 2 + l_23[3]*pow((2*x-1),3))
    def sigma_E31(x):
        return x * (1 - x) * (l_31[0] + l_31[1] * (2 * x - 1) + l_31[2] * (2 * x - 1) ** 2 + l_31[3]*pow((2*x-1),3))
    def sigma_E32(x):
        return x * (1 - x) * (l_23[0] - l_23[1] * (2 * x - 1) + l_23[2] * (2 * x - 1) ** 2 - l_23[3]*pow((2*x-1),3))

    # Calculate eta values
    def integrand1(x):
        return (sigma_E12(x) - sigma_E13(x)) ** 2
    def integrand2(x):
        return (sigma_E21(x) - sigma_E23(x)) ** 2
    def integrand3(x):
        return (sigma_E31(x) - sigma_E32(x)) ** 2

    eta1 = integrate.quad(integrand1, 0, 1)[0]
    eta2 = integrate.quad(integrand2, 0, 1)[0]
    eta3 = integrate.quad(integrand3, 0, 1)[0]

    # Calculate epsilon values
    epsilon12 = eta1 / (eta1 + eta2)
    epsilon23 = eta2 / (eta2 + eta3)
    epsilon31 = eta3 / (eta3 + eta1)
    # print(eta1, eta2, eta3, epsilon12, epsilon23, epsilon31)

    # Calculate X values
    X1 = x1 + x3 * epsilon12
    X2 = x2 + x1 * epsilon23
    X3 = x3 + x2 * epsilon31

    # Calculate w values
    w12 = (x1 * x2) / (X1 * (1 - X1))
    w23 = (x2 * x3) / (X2 * (1 - X2))
    w31 = (x3 * x1) / (X3 * (1 - X3))

    # Calculate G
    G = w12 * sigma_E12(X1) + w23 * sigma_E23(X2) + w31 * sigma_E31(X3)

    # Calculate partial derivatives using finite differences
    h = 1e-10

    def calculate_G(x1_val, x2_val, x3_val):
        # Calculate new X values based on new x values
        X1_new = x1_val + x3_val * epsilon12
        X2_new = x2_val + x1_val * epsilon23
        X3_new = x3_val + x2_val * epsilon31

        w12_new = (x1_val * x2_val) / (X1_new * (1 - X1_new))
        w23_new = (x2_val * x3_val) / (X2_new * (1 - X2_new))
        w31_new = (x3_val * x1_val) / (X3_new * (1 - X3_new))

        return w12_new * sigma_E12(X1_new) + w23_new * sigma_E23(X2_new) + w31_new * sigma_E31(X3_new)

    # Calculate partial derivatives
    dG1 = (calculate_G(x1 + h, x2, x3) - calculate_G(x1 - h, x2, x3)) / (2 * h)
    dG2 = (calculate_G(x1, x2 + h, x3) - calculate_G(x1, x2 - h, x3)) / (2 * h)
    dG3 = (calculate_G(x1, x2, x3 + h) - calculate_G(x1, x2, x3 - h)) / (2 * h)
    # # 验证

    # Calculate deltaGB values
    deltaGB1 = G + dG1 - x1 * dG1 - x2 * dG2 - x3 * dG3
    deltaGB2 = G + dG2 - x1 * dG1 - x2 * dG2 - x3 * dG3
    deltaGB3 = G + dG3 - x1 * dG1 - x2 * dG2 - x3 * dG3

    # print(deltaGB1, deltaGB2, deltaGB3)

    # 北科博士段论文写法，去除第一组元偏导的项
    GB1 = G -x2*dG2 -x3*dG3
    GB2 = G +(1-x2)*dG2 -x3*dG3
    GB3 = G -x2*dG2 +(1-x3)*dG3

    # print(GB1, GB2, GB3)

    # 以上是计算超额吉布斯能的，以下通过多元butler方程计算表面张力

    R = 8.314
    N0 = 6.02E+23
    # 二元 体系
    L = 1.091
    A = [L * pow(N0,1/3) * pow(M[0]/rou[0], 2/3), L * pow(N0, 1/3) * pow(M[1]/rou[1], 2/3), L * pow(N0, 1/3) * pow(M[2]/rou[2], 2/3)]

    def equations(vars):
        xs1, xs2 = vars
        if xs1 <= 0 or xs1 >= 1 or xs2 <= 0 or xs2 >= 1 or (1 - xs1 - xs2) <= 0:
            return [1e6, 1e6]  # 返回大数以表示无效解
        try:
            term1 = theta[0] + (R * T * log(xs1 / x1) + (beta - 1) * GB1) / A[0]
            term2 = theta[1] + (R * T * log(xs2 / x2) + (beta - 1) * GB2) / A[1]
            term3 = theta[2] + (R * T * log((1 - xs1 - xs2) / x3) + (beta - 1) * GB3) / A[2]

            # 方程组
            eq1 = term1 - term3
            eq2 = term2 - term3
            return [eq1, eq2]
        except (ValueError, ZeroDivisionError):
            return [1e6, 1e6]  # 返回大数以表示无效解

    # initial_guess = [0.2, 0.2]
    # solution = fsolve(equations, initial_guess)
    # xs1_solution, xs2_solution = solution
    #
    # solution = root(equations, initial_guess, method='lm')  # 使用Levenberg-Marquardt算法
    from scipy.optimize import root
    initial_guess = [0.1, 0.1]
    solution = root(equations, initial_guess, method='lm')  # 使用Levenberg-Marquardt算法

    if solution.success:
        xs1_solution, xs2_solution = solution.x
        result = theta[0] + (R * T * log(xs1_solution / x1) + (beta - 1) * GB1) / A[0]
    else:
        print("优化未收敛，请检查初始条件")

    return theta[0] + (R * T * log(xs1_solution / x1) + (beta - 1) * GB1) / A[0]
# # #
# def calculate_butler_3_models(modle,T,beta,l_12,l_23,l_13,theta,M,rou,x1,x2,x3):
#     calculator = PE_cal.PE_calculator()  # 创建实例
#     def calculate_G(modle,l_12,l_23,l_13,x1,x2,x3):
#         if modle == "Muggianu":
#             G123 = calculator.Muggianu(x1,x2,x3,l_12[0],l_12[1],l_12[2],l_12[3],l_23[0],
#                                       l_23[1],l_23[2],l_23[3],l_13[0],l_13[1],l_13[2],l_13[3])

#         if modle=="Toop":
#             G123 = calculator.Toop(x1,x2,x3,l_12[0],l_12[1],l_12[2],l_12[3],l_23[0],
#                                       l_23[1],l_23[2],l_23[3],l_13[0],l_13[1],l_13[2],l_13[3])

#         if modle=="Hillert":
#             G123 = calculator.Hillert(x1,x2,x3,l_12[0],l_12[1],l_12[2],l_12[3],l_23[0],
#                                       l_23[1],l_23[2],l_23[3],l_13[0],l_13[1],l_13[2],l_13[3])

#         if modle=="Kohler":
#             G123 = calculator.Koler(x1,x2,x3,l_12[0],l_12[1],l_12[2],l_12[3],l_23[0],
#                                       l_23[1],l_23[2],l_23[3],l_13[0],l_13[1],l_13[2],l_13[3])

#         if modle=="Chou":
#             G123 = calculator.Chou(x1,x2,x3,l_12[0],l_12[1],l_12[2],l_12[3],l_23[0],
#                                       l_23[1],l_23[2],l_23[3],l_13[0],l_13[1],l_13[2],l_13[3])
#         return G123

#     G = calculate_G(modle,l_12, l_23, l_13, x1, x2, x3)

#     h = 1e-12
#     # 偏微分计算
#     dG1 = (calculate_G(modle,l_12,l_23,l_13,x1 + h, x2, x3,) - calculate_G(modle,l_12,l_23,l_13,x1 - h, x2, x3)) / (2 * h)
#     dG2 = (calculate_G(modle,l_12,l_23,l_13,x1, x2 + h, x3) - calculate_G(modle,l_12,l_23,l_13,x1, x2 - h, x3)) / (2 * h)
#     dG3 = (calculate_G(modle,l_12,l_23,l_13,x1, x2, x3 + h) - calculate_G(modle,l_12,l_23,l_13,x1, x2, x3 - h)) / (2 * h)

#     # Calculate deltaGB values
#     deltaGB1 = G + dG1 - x1 * dG1 - x2 * dG2 - x3 * dG3
#     deltaGB2 = G + dG2 - x1 * dG1 - x2 * dG2 - x3 * dG3
#     deltaGB3 = G + dG3 - x1 * dG1 - x2 * dG2 - x3 * dG3

#     # print(deltaGB1, deltaGB2, deltaGB3)

#     # 北科博士段论文写法，去除第一组元偏导的项
#     GB1 = G -x2*dG2 -x3*dG3
#     GB2 = G +(1-x2)*dG2 -x3*dG3
#     GB3 = G -x2*dG2 +(1-x3)*dG3

#     # print(GB1, GB2, GB3)

#     # 以上是计算超额吉布斯能的，以下通过多元butler方程计算表面张力

#     R = 8.314
#     N0 = 6.02E+23
#     # 二元 体系
#     L = 1.091
#     A = [L * pow(N0,1/3) * pow(M[0]/rou[0], 2/3), L * pow(N0, 1/3) * pow(M[1]/rou[1], 2/3), L * pow(N0, 1/3) * pow(M[2]/rou[2], 2/3)]

#     def equations(vars):
#         xs1, xs2 = vars
#         if xs1 <= 0 or xs1 >= 1 or xs2 <= 0 or xs2 >= 1 or (1 - xs1 - xs2) <= 0:
#             return [1e6, 1e6]  # 返回大数以表示无效解
#         try:
#             term1 = theta[0] + (R * T * log(xs1 / x1) + (beta - 1) * GB1) / A[0]
#             term2 = theta[1] + (R * T * log(xs2 / x2) + (beta - 1) * GB2) / A[1]
#             term3 = theta[2] + (R * T * log((1 - xs1 - xs2) / x3) + (beta - 1) * GB3) / A[2]

#             # 方程组
#             eq1 = term1 - term3
#             eq2 = term2 - term3
#             return [eq1, eq2]
#         except (ValueError, ZeroDivisionError):
#             return [1e6, 1e6]  # 返回大数以表示无效解

#     # initial_guess = [0.2, 0.2]
#     # solution = fsolve(equations, initial_guess)
#     # xs1_solution, xs2_solution = solution
#     #
#     # solution = root(equations, initial_guess, method='lm')  # 使用Levenberg-Marquardt算法
#     from scipy.optimize import root
#     initial_guess = [0.33, 0.33]
#     solution = root(equations, initial_guess, method='lm')  # 使用Levenberg-Marquardt算法

#     if solution.success:
#         xs1_solution, xs2_solution = solution.x
#         result = theta[0] + (R * T * log(xs1_solution / x1) + (beta - 1) * GB1) / A[0]
#     else:
#         print("优化未收敛，请检查初始条件")

#     return theta[0] + (R * T * log(xs1_solution / x1) + (beta - 1) * GB1) / A[0]

# # python R-K.py "C:/Users/wang4/OneDrive/热力学软件/Mg-Al-Zn-RE-1873/二元/Zn-Ce-1873K.csv"
# # print(butler(1873,0.75,4044.118,-3729,	2000,0,[0.33,0.7253],[0.024305,0.026982],[1740,2700],0.94))
# #
# # # 生成99个点
# # 生成三元全浓度点并计算表面张力
# def binary_calculatediffT(T,beta):
#     sys = ['Mg_Al', 'Mg_La', 'Mg_Ce', 'Mg_Gd', 'Mg_Y', 'Mg_Nd', 'Al_La', 'Al_Ce', 'Al_Gd', 'Al_Y', 'Al_Nd']
#     x_values = np.linspace(1573, 1873, 100)
#     x2_values = 1-x_values
#     st_values = []
#     all_results = pd.DataFrame()
#     for sys in sys:
#         system = sys
#         L12, theta, M, rou = parameters(system, T)
#         for x in x_values:
#             st = 1000*butler(T,beta,L12,theta,M,rou,x)
#             st_values.append(st)
#             all_results[sys+f''] = x_values
#     df = pd.DataFrame({
#         'x1'+sys: x_values,
#         'x2': x2_values,
#         'st': st_values
#     })
#     return df

# def binary_calculate(T,beta):
#     sys = ['Mg_Al', 'Mg_La', 'Mg_Ce', 'Mg_Gd', 'Mg_Y', 'Mg_Nd', 'Al_La', 'Al_Ce', 'Al_Gd', 'Al_Y', 'Al_Nd']
#     x_values = np.arange(0.001,0.999,0.002)
#     all_results = pd.DataFrame()
#     all_results['x'] = x_values
#     for sys in sys:
#         system = sys
#         L12, theta, M, rou = parameters(system, T)
#         st_values = []
#         exst_values = []
#         for x in x_values:
#             st = 1000*butler(T,beta,L12,theta,M,rou,x)
#             st_values.append(st)
#             exst = st-1000*x*theta[0]-1000*(1-x)*theta[1]
#             exst_values.append(exst)
#         all_results[sys+'-st'] = st_values
#         all_results[sys+'-exst'] = exst_values

#     all_results.to_csv("C:/Users/wang4/OneDrive/热力学软件/Mg-Al-Zn-RE-1873/Mg-Al-RE二元_"+str(T)+"K.csv")


# def binary_calculate_diffT(beta):
#     sys = ['Mg_La', 'Mg_Ce', 'Mg_Gd', 'Mg_Y', 'Mg_Nd', 'Al_La', 'Al_Ce', 'Al_Gd', 'Al_Y', 'Al_Nd']
#     Tem = np.arange(1573, 1873, 1)
#     x_values = [0.01,0.03,0.05] # RE
#     all_results = pd.DataFrame()
#     all_results['T'] = Tem
#     for sys in sys:
#         system = sys
#         for x in x_values:
#             st_values = []
#             for T in Tem:
#                 L12, theta, M, rou = parameters(system, T)
#                 x1 = 1 - x
#                 st = 1000*butler(T,beta,L12,theta,M,rou,x1)
#                 st_values.append(st)
#             all_results[sys+f'x{x:.2f}'] = st_values
#     all_results.to_csv("C:/Users/wang4/OneDrive/热力学软件/Mg-Al-Zn-RE-1873/Mg-RE&Al-RE_1573K-1873K(RE%=1,3,5)).csv")

# # binary_calculate_diffT(0.75)

# def ternary_calculate_all(T, beta):
#     #     # 狄利克雷分布
#     #     points = np.random.dirichlet(alpha=[1, 1, 1], size=1000)
#     #     # 返回结果，每行是一个点[x, y, z]
#     # 生成点
#     compositions = np.random.dirichlet(alpha=[1, 1, 1], size=2000)
#     all_results = pd.DataFrame()
#     x1_values = []
#     x2_values = []
#     x3_values = []
#     st_La = []
#     st_Ce = []
#     st_Gd = []
#     st_Y = []
#     st_Nd = []
#     # 遍历每一个三元成分点，计算对应的st值
#     for comp in compositions:
#         x1, x2, x3 = comp
#         if 0.001 < x1 < 1-1e-3 and 1e-3 < x2 < 1-1e-3 and 1e-3 < x3 < 1-1e-3:
#             x1_values.append(x1)
#             x2_values.append(x2)
#             x3_values.append(x3)
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_La', T)
#             L13, theta3, M3, r3 = parameters('Mg_La', T)
#             st1 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                             [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#             st_La.append(st1)
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_Ce', T)
#             L13, theta3, M3, r3 = parameters('Mg_Ce', T)
#             st2 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                             [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#             st_Ce.append(st2)
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_Gd', T)
#             L13, theta3, M3, r3 = parameters('Mg_Gd', T)
#             st3 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                             [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#             st_Gd.append(st3)
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_Y', T)
#             L13, theta3, M3, r3 = parameters('Mg_Y', T)
#             st4 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                             [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#             st_Y.append(st4)
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_Nd', T)
#             L13, theta3, M3, r3 = parameters('Mg_Nd', T)
#             st5 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                             [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#             st_Nd.append(st5)

#     all_results['Mg'] = x1_values
#     all_results['Al'] = x2_values
#     all_results['La'] = st_La
#     all_results['Ce'] = st_Ce
#     all_results['Gd'] = st_Gd
#     all_results['Y'] = st_Y
#     all_results['Nd'] = st_Nd
#     all_results.to_csv("C:/Users/wang4/OneDrive/热力学软件/Mg-Al-Zn-RE-1873/Mg-Al-RE全浓度"+str(T)+"K.csv", index=False)



# def ternary_calculate_all_modles(modle,T,beta,l_12,l_23,l_13,theta,M,rou):
#     #     # 狄利克雷分布
#     #     points = np.random.dirichlet(alpha=[1, 1, 1], size=1000)
#     #     # 返回结果，每行是一个点[x, y, z]
#     # 生成点
#     compositions = np.random.dirichlet(alpha=[1, 1, 1], size=2000)
#     st_values = []
#     x1_values = []
#     x2_values = []
#     x3_values = []
#     # 遍历每一个三元成分点，计算对应的st值
#     for comp in compositions:
#         x1, x2, x3 = comp

#         st = 1000*calculate_butler_3_models(modle,T,beta,l_12,l_23,l_13,theta,M,rou,x1,x2,x3)
#         st_values.append(st)
            
#     df = pd.DataFrame({
#         'x1': compositions[:,0],
#         'x2': compositions[:,1],
#         'x3': compositions[:,2],
#         'st': st_values
#     })
#     #
#     df.to_csv("C:/Users/wang4/OneDrive/热力学软件/Mg-Al-Zn-RE-1873/"+modle+"-Mg-Al-La-1873K.csv", index=False)



# def ternary_calculate_all_equalstep(T, beta, l_12, l_23, l_13, theta, M, rou):
#     # 生成点
#     x1_values = []
#     x2_values = []
#     x3_values = []
#     st_values = []
#     x = np.linspace(0.001, 0.1, 20)
#     y = np.linspace(0.01, 0.99, 40)
#     # 遍历可能的x,y组合,满足x+y+z=1且所有坐标非负
#     for x1 in x:
#         for x2 in y:
#             x3 = 1-x1-x2
#             if 1e-5 < x3 < 1-1e-5:
#                 st = calculate_butler_3(T, beta, l_12, l_23, l_13, theta, M, rou, x1, x2, x3)
#                 st_values.append(st)
#                 x1_values.append(x1)
#                 x2_values.append(x2)
#                 x3_values.append(x3)
#     for x2 in x:
#         for x3 in y:
#             x1 = 1-x2-x3
#             if 1e-5 < x1 < 1-1e-5:
#                 st = calculate_butler_3(T, beta, l_12, l_23, l_13, theta, M, rou, x1, x2, x3)
#                 st_values.append(st)
#                 x1_values.append(x1)
#                 x2_values.append(x2)
#                 x3_values.append(x3)
#     for x3 in x:
#         for x1 in y:
#             x2 = 1-x1-x3
#             if 1e-5 < x2 < 1-1e-5:
#                 st = calculate_butler_3(T, beta, l_12, l_23, l_13, theta, M, rou, x1, x2, x3)
#                 st_values.append(st)
#                 x1_values.append(x1)
#                 x2_values.append(x2)
#                 x3_values.append(x3)

#     df = pd.DataFrame({
#         'x1': x1_values,
#         'x2': x2_values,
#         'x3': x3_values,
#         'st': st_values
#     })
#     df.to_csv("C:/Users/wang4/OneDrive/热力学软件/Mg-Al-Zn-RE-1873/Mg-Al-Nd-1873K.csv", index=False)



# # 定比例
# def ternary_fixedratio_calculation():
#     x_values = np.linspace(0.01, 0.99, 100)
#     y_values = []
#     for x in x_values:
#         # a = 组元1 ：组元2/3
#         a = 9 / 1
#         y = calculate_butler_3(1873,0.75,
#                                [4044.118,-3729,	2000, 0],
#                                [-60871.19,	-28030,	-6708,0],
#                                [-6248.0468,	-14046.0441,	-7352.9,0],
#                                [0.33,0.7253,0.636024],
#                                [0.024305,0.026982,0.14424],
#                                [1740,2700,7010],
#                                x,(1-x)*(a/(1+a)),(1-x)*(1/(1+a)))
#         y_values.append(y)

#     df = pd.DataFrame(columns=['X','st'])
#     df['X'] = x_values
#     df['st'] = y_values
#     a = str(round(a,2))

#     df.to_csv("C:/Users/wang4/OneDrive/热力学软件/Mg-Al-Zn-RE-1873/三元等比例/Mg-Al-Gd-(Mg-Nd="+a+")-1873K.csv", index=False)

# def ternary_fixedmolefraction_calculate_models(x2,models):
#     x_values = np.linspace(0.001, 0.05, 200)
#     st_values = []

#     for x in x_values:
#         st = 1000*calculate_butler_3_models(models,1873,0.75,
#                                [4044.118,-3729,2000, 0],
#                                [-89252,0,0,0],
#                                [-8775.4435,-3108.1503,	-15106.9],
#                                [0.33,0.7253,0.644912],
#                                [0.024305,0.026982,0.140116],
#                                [1740,2700,6770],
#                                1-x2-x,x2,x)
#         st_values.append(st)
#     return x_values, st_values

# def system3(system3):
#     sys = ['Mg_Al_La', 'Mg_Al_Ce', 'Mg_Al_Gd', 'Mg_Al_Y', 'Mg_Al_Nd']
#     if sys[0] == system3:
#         sys1 = 'Mg_Al'
#         sys2 = 'Al_La'
#         sys3 = 'Mg_La'
#     if system3 == sys[1]:
#         sys1,sys2,sys3 ='Mg_Al','Al_Ce','Mg_Ce'
#     if system3 == sys[2]:
#         sys1,sys2,sys3 ='Mg_Al','Al_Gd','Mg_Gd'
#     if system3 == sys[3]:
#         sys1,sys2,sys3 ='Mg_Al','Al_Y','Mg_Y'
#     if system3 == sys[4]:
#         sys1,sys2,sys3 ='Mg_Al','Al_Nd','Mg_Nd'

#     return sys1,sys2,sys3
# def ternary_fixedmolefraction_calculate(beta,T):
#     Al = [0.02,0.06,0.1]
#     x_values = np.arange(0.001, 0.05, 0.001) # RE
#     all_results = pd.DataFrame()
#     all_results['RE'] = x_values
#     for x2 in Al:
#         st_La = []
#         st_Ce = []
#         st_Gd = []
#         st_Y = []
#         st_Nd = []
#         for x3 in x_values:
#             x1 = 1-x2-x3
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_La', T)
#             L13, theta3, M3, r3 = parameters('Mg_La', T)
#             st1 = 1000*calculate_butler_3(T,beta,L12,L23,L13,[theta1[0],theta1[1],theta2[1]],[M1[0],M1[1],M1[1]],[r1[0],r1[1],r2[1]],x1,x2,x3)
#             st_La.append(st1)
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_Ce', T)
#             L13, theta3, M3, r3 = parameters('Mg_Ce', T)
#             st2 = 1000*calculate_butler_3(T,beta,L12,L23,L13,[theta1[0],theta1[1],theta2[1]],[M1[0],M1[1],M1[1]],[r1[0],r1[1],r2[1]],x1,x2,x3)
#             st_Ce.append(st2)
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_Gd', T)
#             L13, theta3, M3, r3 = parameters('Mg_Gd', T)
#             st3 = 1000*calculate_butler_3(T,beta,L12,L23,L13,[theta1[0],theta1[1],theta2[1]],[M1[0],M1[1],M1[1]],[r1[0],r1[1],r2[1]],x1,x2,x3)
#             st_Gd.append(st3)
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_Y', T)
#             L13, theta3, M3, r3 = parameters('Mg_Y', T)
#             st4 = 1000*calculate_butler_3(T,beta,L12,L23,L13,[theta1[0],theta1[1],theta2[1]],[M1[0],M1[1],M1[1]],[r1[0],r1[1],r2[1]],x1,x2,x3)
#             st_Y.append(st4)
#             L12, theta1, M1, r1 = parameters('Mg_Al', T)
#             L23, theta2, M2, r2 = parameters('Al_Nd', T)
#             L13, theta3, M3, r3 = parameters('Mg_Nd', T)
#             st5 = 1000*calculate_butler_3(T,beta,L12,L23,L13,[theta1[0],theta1[1],theta2[1]],[M1[0],M1[1],M1[1]],[r1[0],r1[1],r2[1]],x1,x2,x3)
#             st_Nd.append(st5)

#         all_results['Mg-Al-'+str(x2)+'La'] = st_La
#         all_results['Mg-Al-'+str(x2)+'Ce'] = st_Ce
#         all_results['Mg-Al-'+str(x2)+'Gd'] = st_Gd
#         all_results['Mg-Al-'+str(x2)+'Y'] = st_Y
#         all_results['Mg-Al-'+str(x2)+'Nd'] = st_Nd

#     all_results.to_csv("C:/Users/wang4/OneDrive/热力学软件/Mg-Al-Zn-RE-1873/Mg-Al0.02&0.06&0.1-RE-"+str(T)+"K(RE%=0-5)).csv",index=False)



# def fumeijiao(T,beta):
#     points = np.random.uniform(low=1e-5, high=0.2, size=(1000,2))
#     x1_values = []
#     x2_values = []
#     x3_values = []
#     st_La = []
#     st_Ce = []
#     st_Gd = []
#     st_Y = []
#     st_Nd = []
#     all_results = pd.DataFrame()
#     for comp in points:
#         x2, x3 = comp
#         x2_values.append(x2)
#         x3_values.append(x3)
#         x1 = 1-x2-x3
#         x1_values.append(x1)
#         L12, theta1, M1, r1 = parameters('Mg_Al', T)
#         L23, theta2, M2, r2 = parameters('Al_La', T)
#         L13, theta3, M3, r3 = parameters('Mg_La', T)
#         st1 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                         [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#         st_La.append(st1)
#         L12, theta1, M1, r1 = parameters('Mg_Al', T)
#         L23, theta2, M2, r2 = parameters('Al_Ce', T)
#         L13, theta3, M3, r3 = parameters('Mg_Ce', T)
#         st2 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                         [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#         st_Ce.append(st2)
#         L12, theta1, M1, r1 = parameters('Mg_Al', T)
#         L23, theta2, M2, r2 = parameters('Al_Gd', T)
#         L13, theta3, M3, r3 = parameters('Mg_Gd', T)
#         st3 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                         [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#         st_Gd.append(st3)
#         L12, theta1, M1, r1 = parameters('Mg_Al', T)
#         L23, theta2, M2, r2 = parameters('Al_Y', T)
#         L13, theta3, M3, r3 = parameters('Mg_Y', T)
#         st4 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                         [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#         st_Y.append(st4)
#         L12, theta1, M1, r1 = parameters('Mg_Al', T)
#         L23, theta2, M2, r2 = parameters('Al_Nd', T)
#         L13, theta3, M3, r3 = parameters('Mg_Nd', T)
#         st5 = 1000 * calculate_butler_3(T, beta, L12, L23, L13, [theta1[0], theta1[1], theta2[1]],
#                                         [M1[0], M1[1], M1[1]], [r1[0], r1[1], r2[1]], x1, x2, x3)
#         st_Nd.append(st5)

#     all_results['Mg'] = x1_values
#     all_results['Ak'] = x2_values
#     all_results['La'] = st_La
#     all_results['Ce'] = st_Ce
#     all_results['Gd'] = st_Gd
#     all_results['Y'] = st_Y
#     all_results['Nd'] = st_Nd
#     #
#     all_results.to_csv("C:/Users/wang4/OneDrive/热力学软件/Mg-Al-Zn-RE-1873/Mg-Al-RE(Al,RE_0-20%)-1573K.csv", index=False)

